#include "func_trapz.h"
#include <malloc.h>


 int func_trapz(complex (*func)(double x, void * context), void * ctx,
		double a, double b,
		double eps, complex * result, complex ** result_arr,
		struct func_trapz_stats * stats){
	if ( eps >= 1 ){
		return 2;
	}

	// Calculating second derrivative maximum
	double fiiabsmax = -1;
	uint32_t nknot_fii;
	int rv = func_secondderr_max(func, ctx, a, b, 0.01, & fiiabsmax, & nknot_fii);
	if(rv){
		fprintf(stderr, "Failed to calculate second derrivative, status=%i\n", rv);
		return 10 + rv;
	}
	// We shal start iterations from some number of knots. Let's use
	// nknot_fii as a guess: at least, we will make calculations
	// of no more difficulty than when calculating second derrivative
	// maximum.
	complex sum;
	uint32_t N = nknot_fii;
	int nitr=0;//iteration counter for staistics
	complex * int_arr = NULL;// array for piece integrals, if requested
	for(;;){//iterations with different N values
#ifdef DEBUG
		printf("Iteration № %i\n", N);
#endif
		nitr += 1;
		double h = (b-a)/(N-1);
		/* renew piece integral array */
		if(NULL != result_arr){
			free(int_arr);
			int_arr = malloc(sizeof(complex) * (N - 1));
			if(NULL == int_arr){
				return 4;
			}
		}

		uint32_t o;
		sum = 0;//integral reset
		complex prev = 0;//Value at previous knot @ o >= 1
		for(o=0;o<N;o+=1){
			double x = a + ((b - a) * o) / (N-1);
			complex here = (*func)(x,ctx);
			if(o>=1){
				complex dint = (prev + here) / 2 * h;
				sum += dint;
				if(NULL != result_arr){
					int_arr [o-1] = dint;
				}
			}
			//push
			prev = here;
		}
		//estimating absolute error
		double delta = 5.0/6.0 * pow((b-a),3) / pow(N,2) * fiiabsmax;
		//and relative error will be
		double nu = delta / cabs(sum); //`delta/approximate_integral_value`
		// To make `delta/true_integral_value` < 2*mu it's enough to check
		// that `nu/approximate_integral_value`< mu providen mu < 1/2
		// We want to limit the whole with eps.
		// mu = eps/2, eps/2 < 1/2 <=> eps<1
		if ( nu < eps / 2 ){
			//Precision achieved
			break;
		}
		//Precision is not enough
		//Which N would we need, if sum were like here
		double newNf = sqrt(5.0/3.0 * pow(b-a,3) * fiiabsmax / cabs(sum) / eps );
		if (newNf > UINT32_MAX){
			free(int_arr);
			return 1;
		}
		uint32_t newN = (uint32_t)newNf + 1;
		if (newN <= N){
			free(int_arr);
			return 3;
		}
		//N
		N = newN;
	}
	//statistics
	if(NULL != stats){
		stats -> nitr = nitr;
		stats -> nknot = N;
		stats -> fii = fiiabsmax;
	}
	if(NULL != result_arr){
		* result_arr = int_arr; //User will free() this
	}
	* result = sum;
	return 0;//Success
}
